# load libraries ----
# reporting
library(knitr)
# visualization
library(ggplot2)
library(ggthemes)
library(patchwork)
# data wrangling
library(dplyr)
library(tidyr)
library(readr)
library(skimr)
library(janitor)
library(magrittr)
library(lubridate)
library(glue)
library(stringr)
# lake models
library(GLMr)
library(glmtools)
# mapping
library(sf)
library(ggspatial)
library(rnaturalearth)
library(rnaturalearthhires)
library(rnaturalearthdata)
# set other options ----
options(scipen=999)
knitr::opts_chunk$set(
tidy = FALSE,
message = FALSE,
warning = FALSE)33 Teleconnections II
Learning Objectives
After completing this tutorial you should be able to
- describe the concepts of macrosystems and teleconnections and how different ecological processes can interact at local, regional, and global scales.
- set up and run lake models to simulate lake temperatures, thermal structure, and ice cover in multiple lakes.
- formulate hypotheses about the effects of teleconnected climate scenarios on different lake modules and test them using simulations.
- describe how local characteristics modify global-scale climate forcing effects on lake temperatures and ice cover.
Download the directory for this project here, make sure the directory is unzipped and move it to your bi328 directory. You can open the Rproj for this module either by double clicking on it which will launch Rstudio or by opening Rstudio and then using File > Open Project or by clicking on the Rproject icon in the top right of your program window and selecting Open Project.
There should be a file named 33_elnino_II.qmd in that project directory. Use that file to work through this tutorial - you will hand in your rendered (“knitted”) quarto file as your homework assignment. So, first thing in the YAML header, change the author to your name. You will use this quarto document to record your answers. Remember to use comments to annotate your code; at minimum you should have one comment per code set1 you may of course add as many comments as you need to be able to recall what you did. Similarly, take notes in the document as we discuss discussion/reflection questions but make sure that you go back and clean them up for “public consumption”.
1 You should do this whether you are adding code yourself or using code from our manual, even if it isn’t commented in the manual… especially when the code is already included for you, add comments to describe how the function works/what it does as we introduce it during the participatory coding session so you can refer back to it.
Before you get started, let’s make sure to read in all the packages we will need. Double check that everything is installed that you need.
33.1 Run & analyze typical El Nino year scenario
Let’s re-define some of our folders that hold our input files necessary for the simulations we’ve set up to run.
# specify your lake
LakeName <- "Sunapee"
# specify your lake simulation folder for observed data
sim_folder <- file.path("sim_lakes", LakeName)
# set filepath for output of your baseline data for your lake
baseline <- file.path(sim_folder, "output.nc")
# specify new sim folder for El Nino year
sim_folder_typical <- glue("sim_lakes/{LakeName}_typical")
# create file pat for extreme el nino simulation
sim_folder_extreme <- "sim_lakes/Sunapee_extreme"Okay, now we’re all set to run our typical El Nino scenario:
run_glm(sim_folder_typical, verbose = TRUE)[1] 0
Before we can look at all the output, we will need to create a file path for our main output file (output.nc) per usual to make it easy to pull all the information we need and create visualizations.
typical_elnino <- file.path(sim_folder_typical, "output.nc")Let’s pull out our temperature for the lake surface and lake bottom and add that to our baseline data. We will want to repeat the same process as above, where we pull up the automatically generated column names2 and then rename them after we join the two data sets3.
2 you will need to check the column names specific to your lake and change the code accordingly!
3 Dates are notoriously difficult to deal with, frequently they might look the same to you, but under the hood there are differences at some level. The code line mutate(DateTime = as.POSIXct(strptime(DateTime, "%Y-%m-%d %H:%M:%S", tz="EST"))) is there to make sure that both DateTime columns are the same class/data type and we can join the data tables
# extract lake depth
LakeDepth <- get_surface_height(baseline)
# get lake temperature for your el nino scenario
LakeTemp_scenario <- get_temp(file = typical_elnino, reference = "surface",
z_out = c(0, min(LakeDepth$surface_height))) %>%
mutate(DateTime = as.POSIXct(strptime(DateTime, "%Y-%m-%d", tz="EST")))
# check column names to be able to rename them downstream
colnames(LakeTemp_scenario)[1] "DateTime" "temp_0" "temp_33.3833669399235"
# combine Lake Temperatures with baseline data
LakeTemp <- get_temp(file = baseline, reference = "surface",
z_out = c(0, min(LakeDepth$surface_height))) %>%
# remember that you will need to use your specific column names to rename your baseline bottom temperature
rename(Baseline_SurfaceTemp = temp_0,
Baseline_BottomTemp = temp_33.3833669399235) %>%
# change date formate so that they match across data frames
mutate(DateTime = as.POSIXct(strptime(DateTime, "%Y-%m-%d", tz="EST"))) %>%
# join with typical el nino year scenario
left_join(LakeTemp_scenario, by = "DateTime") %>%
# rename the column names according to your data set
rename(TypicalElNino_SurfaceTemp = temp_0,
TypicalElNino_BottomTemp = temp_33.3833669399235)Next, let’s pull the ice thickness data for our typical El Nino scenario and add it to the existing LakeIce data.frame.
LakeIce_scenario <- get_var(typical_elnino, "hice") %>%
rename(ice_TypicalElNino = hice) %>%
mutate(DateTime = as.POSIXct(strptime(DateTime, "%Y-%m-%d")))
LakeIce <- get_var(baseline, "hice") %>%
rename(ice_baseline = hice) %>%
mutate(DateTime = as.POSIXct(strptime(DateTime, "%Y-%m-%d"))) %>%
left_join(LakeIce_scenario)Finally, let’s take a look at our lake thermal structure.
plot_temp(file = typical_elnino, fig_path = FALSE)33.2 Run & analyze extreme El Nino scenario
We are ready to run our simulation and create a vector pointing to the output file.
# run simulation
run_glm(sim_folder_extreme)[1] 0
# specify output file
extreme_elnino <- file.path(sim_folder_extreme, "output.nc")Our next step will be to pull out our temperature for the lake surface and lake bottom and add that to our existing LakeTemp data.frame repeating the same process we used above.
LakeTemp_scenario <- get_temp(file = extreme_elnino, reference = "surface",
z_out = c(0, min(LakeDepth$surface_height))) %>%
mutate(DateTime = as.POSIXct(strptime(DateTime, "%Y-%m-%d", tz="EST")))
colnames(LakeTemp_scenario)[1] "DateTime" "temp_0" "temp_33.3833669399235"
LakeTemp <- LakeTemp %>%
mutate(DateTime = as.POSIXct(strptime(DateTime, "%Y-%m-%d", tz="EST"))) %>%
left_join(LakeTemp_scenario, by = "DateTime") %>%
dplyr::rename(ExtremeElNino_SurfaceTemp = temp_0,
ExtremeElNino_BottomTemp = temp_33.3833669399235)Next, let’s pull the ice thickness data for our extreme El Nino scenario and add it to the existing LakeIce data.frame.
LakeIce_scenario <- get_var(extreme_elnino, "hice") %>%
rename(ice_ExtremeElNino = hice) %>%
mutate(DateTime = as.POSIXct(strptime(DateTime, "%Y-%m-%d")))
LakeIce <- LakeIce %>%
mutate(DateTime = as.POSIXct(strptime(DateTime, "%Y-%m-%d"))) %>%
left_join(LakeIce_scenario)Here is an example of where transforming our data can be helpful to better identify patterns. The model gives us the output of ice thickness. However, in many biological and ecological contexts the more helpful metric would be to now what the ice cover is, i.e. how many days of the year is the lake covered in ice.
We can create a new data frame (IceCover) that contains the number of days with ice for each of your scenarios.
IceCover <- LakeIce %>%
pivot_longer(names_to = "scenario", values_to = "ice_thickness", 2:4) %>%
group_by(scenario) %>%
summarize(ice_cover = sum(ice_thickness > 0))A good way to visualize this, would be to plot the relationship of the temperature offset and ice cover duration. This means we need to add the offset to our data.frame. We can do this using a conditional mutate and the helpful function case_when().
typical_offset <- readRDS(glue("results/{LakeName}_typical-offset"))
max_offset <- readRDS(glue("results/{LakeName}_max-offset"))
IceCover <- IceCover %>%
mutate(TempOffset = case_when(scenario == "ice_baseline" ~ 0,
scenario == "ice_TypicalElNino" ~ typical_offset,
scenario == "ice_ExtremeElNino" ~ max_offset))Now it’s straightforward to create plot summarizing the change in ice cover based on the expected temperature offset for a typical and strong El Nino year.
ggplot(IceCover, aes(x = TempOffset, y = ice_cover)) +
geom_line() +
geom_point(aes(fill = scenario), shape = 21, size = 3) +
scale_fill_manual(values = c("darkblue", "red", "darkorange")) +
labs(x = "Temperature offset [C]", y = "ice cover duration [days]")Finally, let’s take a look at our lake thermal structure across our three scenarios.
We’ve previously relied on the glmtools::plot_temp() function. However, down the line you might want to plot a heat plot for another data set down the line. Let’s learn how to use ggplot to plot a heatmap.
First, we need to pull the temperature data for all depths. If we don’t specify z_out as we have done previously to get just the surface and/or bottom temperatures the get_temp() function will return data for all depths.
therm_struct <- get_temp(file = extreme_elnino, reference = "surface") %>%
mutate(DateTime = as.POSIXct(strptime(DateTime, "%Y-%m-%d", tz="EST")))To plot a heatmap we will us the function geom_tile(). We need a data set with a column for date which we will plot on the x-axis (this already exists as DateTime), a column with depth which we will plot on the y-axis, and a column with temperature (this is the one we will color code the tiles by).
This means we need to create a tidy data set. Currently, the column names for each depth are in the format temp_depth. We can use the function stringr::str_sub()4. Your data.frame might have a different number of columns than our example. We can use ncol(therm_struct) to get the number of columns in the data set and use that to specify which columns to gather.
4 stringr is a package in the tidyverse designed to manipulate strings.
therm_struct <- therm_struct %>%
pivot_longer(names_to = "depth", values_to = "temperature", 2:ncol(therm_struct)) %>%
mutate(depth = str_sub(depth, 6), # remove first five characters
depth = as.numeric(depth)) %>% # convert to numeric
filter(!is.na(temperature)) # remove missing values.Now we can plot our heatmap showing the thermal structure for our extreme El Nino year.
We are going to flip the y-axis using scale_y_reverse() so that our surface temperature is at the top of the plot. We can specify the range of values using the limits and breaks options in the scale_fill_viridis_c() function that we are using to specify the color scale. The temperature cannot drop below 0 but your should check the maximum value for your temperature to determine the highest value for your scale5.
5 in my case it is around 27C so I chose 30
ggplot(therm_struct, aes(x = DateTime, y = depth, fill = temperature)) +
geom_tile() +
scale_fill_viridis_c(option = "plasma", limits=c(0, 30), breaks=seq(0, 30, by = 5)) +
scale_y_reverse() +
labs(x = "date", y = "depth [m]")33.3 Acknowledgments
These activities are based on the EDDIE Teleconnections module.6
6 Farrell, K.J. and C.C. Carey. 18 May 2018. Macrosystems EDDIE: Teleconnections. Macrosystems EDDIE Module 3, Version 1.